source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/ADRB2-FOH/'

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg')

sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj$orig.ident |> unique()

# All Adr genes -------
adr.mm.gene <- sobj |> rownames() |> str_subset('^Adr[a-c]') |>
  str_sort()

sobj |> filter(orig.ident == 'PBS') |>
  DotPlot(adr.mm.gene, cols = 'RdYlBu', group.by = 'manual_main') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Adrenoceptors in mouse skin')

adr.mskin.expr <- last_plot() |>
  pluck('data') |>
  as_tibble()

g1 <- adr.mskin.expr |>
  tidyplot(x = features.plot, y = id, color = avg.exp.scaled) +
  geom_point(aes(size = pct.exp))

g2 <- g1 + theme_bw() 

adr.mskin.expr |>
  BubblePlot() +
  labs(title = 'Adrenoceptors in mouse skin') +
  theme_bw(base_size = 6) +
  RotatedAxis() +
  theme(legend.key.size = unit(2, 'mm'), legend.position = 'bottom')

publish_source_plot('mice.skin.adr.genes.dotplot')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot('Adrb2', cols = 'RdYlBu', group.by = 'manual_main')

sobj |>
  filter(orig.ident == 'PBS') |>
  mutate(manual_main = fct_relevel(manual_main, 'Keratinocytes')) |>
  bill.violin('Adrb2', group.by = manual_main) +
  RotatedAxis() +
  NoLegend() +
  labs(title = 'Healthy mice skin Adrb2 expression',
       x = 'Cell type', y = 'Normalized expression')

## Dotplot in KC ---------
sobj.kc <-
  read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj.kc |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c('Adrb2', late.kc, 'Krt5','Krt14'), group.by = 'seurat_clusters',
          cluster.idents = T, cols = 'RdYlBu') +
  labs(title = 'Adrb3 & KC differentiation markers in mice KC',
       x = 'Genes', y = 'KC clusters')

mkc.adrb2.diffmark <- last_plot() |>
  pluck('data')

# TODO: rewrite in tidyplots
mkc.adrb2.diffmark |>
  filter(features.plot == 'Adrb2') |>
  mutate(order = fct_reorder(id, avg.exp), id, .keep = 'none') |>
  right_join(mkc.adrb2.diffmark) |>
  mutate(id = order) |>
  BubblePlot(size = 5) +
  labs(title = 'Adrb2 & KC differentiation marker in mice KC',
       y = 'KC clusters')

sobj |>
  DotPlot(list(Adrb2 = 'Adrb2',
               Late.KC = late.kc,
               Early.KC = c('Krt5','Krt14')), cols = 'RdYlBu')

sobj.kc |>
  DotPlot(c('Adrb2', map(cc.genes.updated.2019, str_to_title)))


## define Adrb2-hi KC ---------
sobj |>
  get_abundance_sc_wide('Adrb2') |>
  filter(Adrb2 > 0) |>
  ggplot(aes(Adrb2)) +
  geom_density()

mkc.adrb2.diffmark |>
  filter(features.plot == 'Adrb2') |>
  ggplot(aes(avg.exp, avg.exp, label = id)) +
  geom_point() + geom_text_repel(nudge_y = .2)

## Adrb2-hi vs lo KC DEGA -------
mb2h.kc.deg <- sobj.kc |>
  FindMarkersAcrossVar(split.by = 'orig.ident', group.by = 'seurat_clusters',
                       ident.1 = c(1,11,3))

mb2h.kc.deg |>
  write_csv('mission/ADRB2-FOH/results/mice.b2hvl.kc.deg.csv')

mb2h.kc.deg |>
  filter(cluster == 'PBS', !str_detect(gene, 'Rik$|^Gm|^AC')) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  plot_bill_volc(group1 = 'Adrb2-high KC',
                 group2 = 'Adrb2-low KC')

library(clusterProfiler)
pbs.mb2h.gsego <- mb2h.kc.deg |>
  filter(cluster == 'PBS', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db', keyType = 'SYMBOL')

pbs.mb2h.gsego <- pbs.mb2h.gsego |>
  simplify()

pbs.mb2h.gsego@result |>
  write_csv('mission/ADRB2-FOH/results/mice.b2h.kc.gsego.csv')

pbs.mb2h.gsego@result |>
  filter(NES > 0, ONTOLOGY == 'BP') |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES',
       title = 'Upregulated GO pathways\nAdrb2-high KC vs Adrb2-low KC')

pbs.mb2h.gsego@result |>
  filter(NES < 0, ONTOLOGY == 'BP') |>
  plot_enrichment(metric = NES, base_col = 'blue') +
  labs(x = 'NES',
       title = 'Downregulated GO pathways\nAdrb2-high KC vs Adrb2-low KC')
